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Abstract The frequency-domain fast boundary element 
method (BEM) combined with the exponential window 
technique leads to an efficient yet simple method for 
elastodynamic analysis. In this paper, the efficiency of 
this method is further enhanced by three strategies. 
Firstly, we propose to use exponential window with 
large damping parameter to improve the conditioning 
of the BEM matrices. Secondly, the frequency domain 
windowing technique is introduced to alleviate the se- 
vere Gibbs oscillations in time-domain responses caused 
by large damping parameters. Thirdly, a solution ex- 
trapolation scheme is applied to obtain better initial 
guesses for solving the sequential linear systems in the 
frequency domain. Numerical results of three typical 
examples with the problem size up to 0.7 million un- 
knowns clearly show that the first and third strategies 
can significantly reduce the computational time. The 
second strategy can effectively eliminate the Gibbs os- 
cillations and result in accurate time-domain responses. 
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1 Introduction 

Transient analysis is encountered in almost every field 
in science and engineering. For most real-world prob- 
lems, numerical simulation is the only viable approach 
for obtaining accurate solutions. The boundary element 
method (BEM) has been well-recognized as a powerful 
numerical method to treat wave propagation problems. 
Its reduced dimensionality, high accuracy and its abil- 
ity in capturing rapid transitions and steep gradients 
of the fields have made this method particularly suit- 
able for problems with complex geometry, for example, 
dynamic analysis of porous solids and fracture anal- 
ysis. In view of the increasingly complex systems that 
we have encountered nowadays, the development of effi- 
cient BEM for large-scale wave propagation simulations 
has become indispensable. 

According to the different approximation solution 
strategies in time space, the BEM for treating elasto- 
dynamic wave propagation problems generally follows 
two approaches, namely, time-domain approaches and 
frequency-domain approaches (see e.g. the reviews by 
Beskos fT] and Costabel 0). Time-domain approaches 
can be further classified into time-stepping methods 
and the space-time integral equation method. In these 
methods the physical problems are directly solved in 
the real time domain, thus one can observe the phe- 
nomenon as it evolves. However, such methods require 
an adequate choice of the time step size. An improper 
time step could lead to instability or numerical damp- 
ing. For recent development of time-domain methods, 
see e.g. HHE]. 

The frequency-domain approach based on the fast 
Fourier transform offers an attractive method for tran- 
sient analysis due to its stability and simplicity [5EIB] ■ 
In this approach, only the frequency-domain equation 
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needs to be solved, which by itself has a broad range 
of applications such as harmonic wave propagation in 
solids pQ, air damping on micro resonators |10j . heat 
transfer [9], to name a few. 

There are two issues need to be dealt with when 
computing the transient response using the inverse Fourier 
transform. The first issue is regarding to the discretiza- 
tion and truncation errors introduced by the discretized 
Fourier transform (DFT) [T7]- The discretization in DFT 
causes wrap-aroud (or aliasing) error. For damped sys- 
tems this error can be removed by adding trailing ze- 
roes to damp out the free vibration before the end 
of the period. However, this scheme is inefficient for 
low damped systems, and more importantly, it does 
not work for non-dissipativc systems such as the elas- 
tic solids considered in this work. To circumvent this 
problem, the exponential window method (EWM) [TH1 
HI?] was employed in the BEM transient clastodynamic 
analysis p2Ull21j . This method adds an artificial damp- 
ing to the system to effectively damp out the free vi- 
brations. By properly choosing the damping parame- 
ter, denoted by 77, in the EWM, the response period T 
can be arbitrarily set independent of the actual time- 
domain response. In addition, it has been found that 
a large 77 can improve the conditioning of the BEM 
coefficient matrix. However one cannot arbitrarily in- 
crease 77 because it also amplifies the truncation error 
of DFT and therefore causes unacceptable Gibbs os- 
cillations in the later time response. In this paper, a 
frequency-domain windowing technique is introduced 
to eliminate the Gibbs oscillations, and consequently 
accurate time response can be obtained in all simula- 
tion periods. The frequency-domain windowing tech- 
nique in combine with the EWM leads to a simple but 
efficient modified Fourier transform approach for BEM 
clastodynamic transient analysis. 

The second issue concerns the computational ef- 
ficiency. In the frequency-domain approach the tran- 
sient problem is transformed into N u independent time- 
harmonic problems. The computational cost grows al- 
most linearly with N u . For many practical problems 
N u > 100 sampling frequencies are often required to 
obtain accurate time-domain solutions. This implies a 
huge computational cost because solving one frequency- 
domain eqation in a three-dimension domain using clas- 
sic BEM is already computationally very costly. To 
resolve this issue, a series of work has been recently 
published on the development of efficient integral equa- 
tion approaches for frequency-domain elastodynamics. 
Examples include, but are not limited to, fast multi- 
pole [THI12U13] . "H-matrix |H] and prccorrected-FFT 
(pFFT) accelerated BEM [15], as well as efficient BEM 
with adaptive cross approximation [IB]. In this paper, 



the pFFT technique is employed to accelerate the frequency- 
domain BEM analysis [3T] due to its ease in implemen- 
tation and high computational efficiency. 

A key factor that affects the computational cost of 
the frequency domain approach is the number of itera- 
tions required to solve each linear system. In order to 
achieve further speedup, a least square extrapolation 
method is proposed to obtain better initial guesses for 
the iterative solution procedure. 

The basic setting of the frequency-domain approach 
for BEM transient analysis is briefly recalled in Section 
[21 The modified Fourier transform method for tran- 
sient elastodynamic anaylsis is proposed in Section[3] In 
Section [4] the precorrected-FFT and the extrapolation 
method for initial guess arc presented. The performance 
of the methods put forward in this paper is tested by 
three examples in Section [5J and a summary is drawn 
in Section [B] 

2 Frequency domain BEM for transient analysis 

2.1 Frequency domain boundary element method 

The first step in the frequency-domain approach for 
transient elastodynamic analysis is to transform the 
time domain into the frequency domain via Fourier 
transform. Let fl G R 3 denote the region of space occu- 
pied by a three-dimensional elastic solid with isotropic 
constitutive properties defined by Lame constants A 
and fi, Poisson's ratio v and mass density p. In a mixed 
boundary value problem, the boundary T = df2 con- 
sists of a displacement boundary r u and a traction 
boundary T a such that F = r u U T a . The displacement 
governing equation for linear elasticity in the frequency 
domain is formulated as 

/iAu(x, to) + (A + /i)VV • 7t(x, uj) = puj 2 u(x, uj), x G fi, 

u(x, uj) = u(x,oj), x G r u , 

a(x,oj) = (& x u)(x,u)) = a(x,oj), x G F a , 

(1) 

where, uj is the circular frequency; u — (u±, ^2,1^3) is the 
displacement in the frequency domain; V and A denote 
the Nabla and the Laplace operators, respectively; 27^ is 
the traction operator defining the stress-strain relation 
based on Hooke's law; u and a are prescribed boundary 
displacement and traction on r u and F a , respectively. 

The boundary integral equation corresponding to 
equation ([T]) is formulated as 

Cij(x)uj(x,u)) + (P.V.) / Ty(x,y;w)iij(y,u;)dr y 

, r (2) 

= J U ij (x,y;u)a j (y,u)dr y , xef 
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where, (P.V.) indicates a Cauchy principal value (CPV) 
of the singular integral; the free-term Cij(x) is equal 
to 0.55ij for a smooth boundary at x; C/y ;(x, y; w) and 
Tij(x, y;u>) denote the jth components of the displace- 
ment and traction generated at y by a unit point force 
applied at x along the direction i, given by 



Uij(x,y;uj) 



1 

4nfx 



[<p(r)dij - ipnrj] , 



T 4j (x,y;w) = G ikj (x,y;uj)n k (y), 
where, 



(3) 



(4) 
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n(y) = (hi, ri2, H3) is the unit normal to J 1 at y directed 
outwards of i7: 



l x ~y| 



are the respective wavenumbers of S and P 
elastic waves, so that 



UJ UJ 

, fop 
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The boundary integral equation ^ is numerically 
solved using a surface discretization with N c triangular 
boundary elements. Piecewise-constant approximation 
is employed in this work. A standard collocation ap- 
proach yields a discretized system of the BIE © of the 
form 



Tu = Ut, 



(6) 



where matrices T and U with dimensions of N by N, 
N = 3N C , consist of boundary integrals of kernel func- 
tions Tij and Uij on each element respectively; the N- 
vectors u and t contain the nodal displacement and 



traction at each element. Note that the contribution of 
the free-term Cy (x) has been absorbed into the matrix 
T. By enforcing the boundary conditions, the linear 
system ([6|) becomes 



A(w)a(w) = b(w) 



(7) 



where, the TV-vector a collects the unknown nodal dis- 
placement and traction components and can be ob- 
tained by solving the above linear system. All the ma- 
trix and vectors in ((7]) are functions of frequency ixs. 
The computational cost for an iterative solver is 0(N 2 ). 
This can be greatly reduced by using the pFFT method 
in section |4~T1 



2.2 Conventional frequency-domain approach 

Equation (TT|) or ^ is solved at a series of discrete fre- 
quencies. The resulting displacements and tractions are 
then transformed back to the time domain by using the 
inverse DFT to obtain the time-domain response. 

Let Auj and At be the circular frequency and time 
resolutions, respectively; and T be the time period of 
the transient response. Given the number of sampling 
points N u , one has the basic relations T = 2it/Auj and 
At = T/N u . Note that the Nyquist frequency /Nyq = 
N LU Auj/(4:Tr) should be chosen such that responses cor- 
responding to frequencies higher than /Nyq are insignif- 
icant and can thus be discarded. 

The procedure for obtaining the transient responses 
using the frequency-domain approach can be summa- 
rized as: 

(1) Determine a frequency resolution Auj and the num- 
ber of sampling points N u . Auj needs to be small 
enough to minimize the loss of useful frequency in- 
formation. 

(2) Sample boundary condition P(t), which could be 
the given boundary displacement or the traction, 
and perform DFT to obtain the frequency-domain 
boundary conditions P(ui) at frequencies {u) k = 
kAu, (* = 0,1,..- ,N U -1)}, 



P(w fc ) 



1 



N u -l 

E 

71=0 



P(nAt)e 



(8) 



(3) Conduct analysis for the first (-^=r + 1) frequencies, 
that is, {uj k = kAuj, (k = 0, 1, • ■ ■ ,N u /2)}, to ob- 
tain the frequency-domain response R(uj) for the 
first (-^ + 1) frequencies. The responses R(uj) cor- 
responding to the last (^r — 1) can be determined 
utilizing the conjugate symmetric property of the 
DFT as, 



R(k)=R*{N u -k) 1 k = N u /2+l, 



1, 
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where, R(k) = R(kAuj) and R* denotes the complex 
conjugate of R. 
(4) Perform IDFT for R(uj) to obtain time-domain re- 
sponse R(t), i.e., the desired boundary displacement 
and/or traction, 



R{nAt) = ^K)e 27r 



(9) 



fe=0 



Discretization and truncation errors arise when the 
forward and inverse Fourier transforms are numerically 
evaluated using DFT and IDFT in © and ©. The 
discretization error is also known as aliasing error. To 
eliminate it, the time-domain signal is required to be 
zero for t > T . The truncation error is due to the trun- 
cation of the frequency spectrum. It causes high fre- 
quency oscillations near the discontinuities of the sig- 
nals, known as Gibbs oscillations. 

In next sections, techniques for reducing the two 
errors are introduced. Consequently, a general and ef- 
ficient modified Fourier transform (MFT) method for 
elastodynamic transient analysis is obtained. The MFT 
method itself is not new. It was first developed for an- 
alyzing power system transients by Mullincux's group 
about four decades ago; see [T7J and the references therein. 
In this paper, it is applied in elastodynamic transient 
analysis. 



3 Modified Fourier transform method 

In the MFT method, two data processing techniques, 
one with respect to the time-domain data and another 
to the frequency-domain data, are introduced to reduce 
the discretization and truncation errors associated with 
the DFT and IDFT. 



3.1 Discretization error & time-domain damping 

As aforementioned, the periodic nature of the DFT 
requires that the time-domain response tends to zero 
at the end of a period. This however often cannot be 
satisfied in the dynamic analysis using the frequency 
domain approach. A common way to circumvent this 
problem for a dissipated system is to add trailing zeroes 
to prolong the duration so that the free vibration can 
be damped out before the end of the period. Such an 
approach is clearly not applicable for a non-dissipated 
system. Another more general method suitable for all 
systems is the exponential window method (EWM) [TBI 
119) . In this method, a complex frequency shift, amount- 
ing to an artificial damping, is introduced to efficiently 
damp out the free vibrations within the selected period. 



The actual response is obtained after a proper scaling 
of the damped response. 

Recently, the EWM was used in 2D and 3D BEM 
elastodynamic transient analysis f20ll2~T] . Accurate time 
domain response can be obtained by a proper choice 
of the damping parameter n > 0. To be specific, let 
tt(x, t) be the displacement in the time domain and 
its Fourier transform is tt(x, w) in ([1]). In applying the 
EWM, the displacement u is scaled by an exponential 
window function e~ r,t to obtain a damped displacement 



He 



i.e. 



e 



-</' 



u(x, t). 



(10) 



Since u is a bounded function, the damped response 
ii cw tends to be zero for large t. The Fourier transform 
of M ew is it ew ( x , k>) = u(x,a> — uq), meaning that u cw 
can be obtained by simply replacing the frequency to in 
u by to — it], where i = y— T. To obtain the transient 
response, one needs to solve the problem at a series of 
complex frequencies uJk = kAaj—ir], k = 0, 1, • • • , N u /2. 
Once the damped response u ow is obtained, the actual 
response u can be obtained as 



u(x, t) 



e vt u evf (x,t). 



(11) 



The choice of the damping parameter t] determines 
the performance of the EWM. On one hand, it can- 
not be too large since the factor e 7,t in (flT|) acts as 
an amplifier which magnifies the cumulative errors due 
to the truncation of the frequency space which will be 
discussed in section 13.21 and BEM error. On the other 
hand, given the factor e~ nt is used to damp the time- 
domain response, it could be supposed that a high value 
for t] is required to sufficiently eliminate the aliasing er- 
ror. A compromise between the two ends leading to the 
following rulc-of-thumb for the selection of t] |18 



k In 10 
n=—^ = KAf In 10. 



(12) 



This rule was used in the work described in 821, [217]. It 
was found in [21] that in BEM analysis a choice of the 
coefficient k between 2 < n < 3 could often lead to 
good time-domain results. 

In addition to reducing the aliasing error, another 
important effect of rj lies in that a positive 77 can im- 
prove the conditioning of the system matrices in ([5]) and 
l[7|) . This is a nature result of the Gershgorin circle theo- 
rem, because the exponential function e~ nt reduces the 
absolute values of the non-diagonal entries of the matri- 
ces, thus making the eigenvalues more clustered. Hence, 
the damping factor e~ nt actually acts as a percondi- 
tioncr (called artifical damping perconditioner, ADP, in 
this paper) whose performance can be enhanced by in- 
creasing 77. As such, a larger n is always desirable. 



Efficiency Improvement of the frequency-domain BEM for Rapid transient elastodynamic analysis 



5 



Matrix preconditioning is a crucial problem in BEM 
analysis. This is particularly true for high frequency and 
multi-domain problems, in which the coefficient matri- 
ces tend to be much more ill-conditioned. Various pre- 
conditioners exist (e.g., [221I23] ) and can be used. How- 
ever, it is well known that the performance of many ex- 
isting preconditioners is generally problem-dependent 
and there always is a trade-off between their perfor- 
mance and computational cost. On the contrary, the 
ADP is simple yet efficient, because its implementation 
requires no extra computation and memory. Whereas, 
it can largely increase the convergence speed of the it- 
erative solvers; see Section [5] 

Meanwhile, the ADP can be used in combination 
with other preconditioners. For instance, in this work 
a block diagonal preconditioner is used together with 
the time-domain damping technique described in this 
section to solve ([7]). 



3.2 Truncation error & frequency-domain windowing 



3.3 A summary 

To summarize, the time-domain damping in conjunc- 
tion with the frequency-domain windowing techniques 
leads to a MFT method for elastodynamic transient 
analysis. The basic implementation procedure of the 
MFT method is similar to that described in section 
[2721 The distinction lies in that in the MFT method 
the sampling frequencies are uik = kAu — ltj instead of 
u) k = kAuj, and the DFT formula © and IDFT formula 
© are replaced by 



N u -1 



p ( u ") = aT £ e-^ At P(nAt)e 



■2-n-ink/N^ 



n=0 



and 



R{nAt) = e r ' nAt £ W^R^e 2 * 



(13) 



(14) 



fc=0 



respectively. 



It is observed that choosing k > 2 in BEM analysis can 
cause unacceptable Gibbs oscillations in time-domain 
responses at later times |21j . Here, a frequency-domain 
data windowing technique is introduced to alleviate the 
the Gibbs oscillations. Consequently, a larger value of k 
can be set. This implies a considerable improvement in 
the conditioning of the BEM matrices when compared 
with k = 2; see evidences in Section [5] 

The implementation of this technique is rather sim- 
ple [17] . Before the damped frequency-domain response 
u ew (x, lo) is transformed back into the time domain 
using IDFT, it is first multiplied by a suitable data 
window W(uj); that is, the function H^(aj)M ew (x, w) is 
transformed in lieu of u ew (x, w) itself. The function of 
the frequency-domain windowing is to reduce the Gibbs 
oscillations by averaging out the oscillations within one 
period. 

There exist a number of window functions. As far 
as BEM elastodynamic analysis is concerned, we find 
that the Hanning and Blackman windows often yield 
satisfactory results. The discretized window functions 
are given by 



Hanning: W(uJk) = 0.5 



1 + cos 



27rfc 



(27tfc 
— — 

/4nk\ 
+ 0.08 cos — - , 

fc = 0,l,--- ,N U -1. 



4 Rapid solution of linear systems 

The main computational cost for the frequency-domain 
approach is the solution of linear systems ([7]) for a se- 
quence of sampling frequencies i-e., the parameter- 
ized linear systems 



A fc a fc 



1... -^L 



2 



(15) 



where, = A(ujk), &fc and b& are defined analogously. 
In this section, two methods are described to consider- 
ably reduce the computational cost, namely, a pFFT 
method for accelerating the BEM matrix-vector multi- 
plication, and a frequency extrapolation method for ob- 
taining good initial guesses and thus reducing the total 
number of iterations in solving the sequential systems. 



4.1 PFFT accelerated elastodynamic BEM 

The computational cost of a matrix-vector multiplica- 
tion in traditional BEM is 0(N 2 ), which could be pro- 
hibitively expensive for large-scale elastodynamic simu- 
lations. Within the last two decades, significant progress 
has been made in the development of accelerated BEM 
techniques. With the accelerated techniques, the com- 
putational complexity, including both the CPU time 
and memory usage, can be considerably reduced, there- 
fore enabling BEM-based large-scale computation. 

The first accelerated technique used in elastody- 
namic BEM is the FMM [53] , which is further developed 
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by several groups; see, e.g., [2511131112] . In [T3] the H- 

matrix method is applied in frequency-domain elasto- 
dynamic BEM. Recently the pFFT method was used to 
accelerate elastodynamic BEM analysis by the present 
authors [21]. The pFFT method is very easy to imple- 
ment and to be parallelized, so a high computational 
efficiency can be achieved with ease. For more informa- 
tion about the pFFT technique and the details of its 
implementation, the readers are referred to [21j . 

Theoretically, for problem domains with small surface- 
to- volume ratios such as porous solids, optimal order 
0(N log N) for CPU time and O(N) for memory can 
be achieved by pFFT, which is comparable to the FMM 
and the H-matrix method. For other problems, the CPU 
time of pFFT is between 0(N log N) and ©(iV 15 logiV). 
In these cases, the FMM and the H-matrix method 
would be advantageous. However numerically, due to 
the small overhead, the actually computational cost of 
pFFT could be less than that of FMM and the U- 
matrix method for problems with N up to order O(10 5 ); 
see [21] for comparisons. 

4.2 Solution extrapolation method 

Iteratively solving each linear system in the sequence 
(fT5|) needs an initial guess, denoted by . Typically, 
aj^ is independently chosen as zero or one vectors. Here 
a least square extrapolation method is proposed to de- 
termine a better initial guess a^ for the solution of a*, 
based on the calculated solutions at previous frequen- 
cies in the sequence; this approach is termed as the 
solution extrapolation method (SEM). 

The motivation is inspired by the fact that solution 
a(w) depend on the frequency oj. If a(w) is a smooth 
function of uj and that the frequency step Aoj is suf- 
ficiently small, it should be possible to obtain a^ by 
extrapolation from the previous K solutions; that is, 

K 

a[ 0) = ^2 s * a k-i = Xs, (16) 

i=l 

where, X = (a fe _i, • • • ,a k ^ K ), s = (si, • • • 7 s K ) T ■ Let 
Y = AfcX. Then the coefficient vector s can be ob- 
tained by minimizing the square error 

||A fc a£ 0) -b fc || 2 = ||Ys-b fc || 2 . (17) 

By computing the QR-dccomposition Y = QR, one 
gets the least square solution of ([TTf 

s = R iQbfc. 

Thus, the initial guess is can be computed as 

af =XR- 1 Qb fe . (18) 



The main computational cost for evaluating a), is the 
computation of K matrix-vector products A&X. 

Optimal value for the number of preceding solutions 
K in (fTB"]) clearly depends on the smoothness of the 
solutions a(w) and thus is problem dependent. For the 
cases studied in this work, it is found that K = 3,4 
is optimal. Generally, a large K does not improve the 
accuracy of the initial guess, because the solution a(ui) 
is often not smooth over a large frquency range and the 
frequency step Aw may not be so small as required. 



5 Numerical examples 

In this section, three examples are selected to validate 
and demonstrate the performance of the MFT method 
in Section[3]and the SEM for initial guess in Section l4T2"l 
The first example is a classical benchmark for transient 
elastodynamic analysis: a prismatic rod subject to a 
sudden applied uniform pressure at one end. The ana- 
lytic solution is available and is employed to evaluate 
the accuracy and the efficiency of the method. The sec- 
ond case is an aluminum plate subject to an impact 
loading. This example is chosen to test the method for 
problems with a time-dependent arbitrary loading. Fi- 
nally to test the ability of the method for solving large- 
scale problems, a unit cube containing 64 spherical cav- 
ities subject to a sudden applied uniform pressure is 
simulated. 

In the following simulations, constant elements are 
used. The pFFT method in [21] is used to accelerate the 
BEM analysis. All the linear systems are preconditioned 
using the block diagonal preconditioncr and solved by 
GMRES with convergence tolerance 10 -5 . Calculations 
are performed on a workstation with a Xeon 3.0GHz 
CPU and 30GB RAM. 



5.1 Example 1: prismatic rod 

As the first test case, a prismatic rod with dimensions 
of 3m x lm x lm as sketched in Figure [1] is simulated. 
The rod is fixed at its left end and subject to a Heav- 
iside traction p(t) = PoH(t) at its right end, where 
Pq = — 10 6 N/m 2 and H(t) is the Heaviside function. 
By setting the Poisson's ratio of the material to zero, 
the 3-D problem is reduced to a 1-D problem, of which 
the analytical solution is available, see [21]. The Young's 
modulus E and the density p of the material are set to 
be E = 2.11 ■ 10 11 Pa and p = 7.85 ■ 10 3 kg/m 3 , respec- 
tively. 
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Fig. 1 Prismatic rod subject to a step traction pulse 

5.1.1 Effects of the ADP and the SEM 

The boundary of the prismatic rod is partitioned into 
12800 triangular elements. The response period and the 
number of samplings are set to be T = 0.0155s and 
N u = 256, respectively. To test the effects of the damp- 
ing coefficient k and the SEM on accelerating the con- 
vergence of the linear system solver, simulations arc 
carried out under the three cases specified by the first 
three columns of Table [TJ For cases without using the 
SEM, the initial guesses are set to be zero vectors. 



Table 1 Three test cases and their CPU times 



Cases 


K, 


SEM 


Total its 


Total time (h) 


Case 1 


2 


X 


29237 


16.0 


Case 2 


4 


X 


21758 


14.4 


Case 3 


4 


/, K = 4 


17562 


12.1 



The number of iterations for solving the sequential 
linear systems (fT5")) corresponding to the first 129 sam- 
pling frequencies are plotted in Figure[5] The total num- 
ber of iterations and the total CPU time for solving all 
the 129 systems in each of the three cases are listed in 
the last two columns of Table [TJ One can see that by 
using a larger damping with k = 4, a 25% reduction 
in the total number of iterations can be achieved when 
compared with the case of K = 2. By further using the 
SEM, the reduction is up to about 40%. However, the 
reduction of the total CPU time is not as considerable 
as the number of iterations. This is because the matrix- 
vector multiplication in this problem is relatively fast 
and the preprocessing time (i.e., time for computing the 
near-field intcrations) occupies a large portion. 

5.1.2 Effects of the frequency windowing technique 

As described in Section 13.21 the frequency windowing 
technique can be employed to eliminate the Gibbs os- 
cillations in the time-domain response. Here this is ver- 
ified using the results of case 3, i.e. n = 4, in Table [TJ 
Our numerical experience shows that both the Black- 
man and Hanning windows obtains almost the same 
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results. Thus only the results of Hanning window are 
presented here. 

Figure 3 and 4 demonstrate both the analytical and 
simulated displacements at the free end and the trac- 
tions at the fixed end as functions of time; the non- 
dimensional time (t ■ c/L, with c and L being the wave 
speed and length of the rod, respectively) is employed in 
the plots. It is observed that, before using the frequency 
windowing technique, both the displacement and trac- 
tion responses show drastic oscillations in later times. 
The magnitude of the traction oscillations is about 4 
times larger than the largest traction value; see Fig- 
ure [4] (a). But after the frequency windowing technique 
is used, the red curves obviously show that the oscilla- 
tions arc effectively removed and the resulting responses 
agree very well with the analytical solutions in all the 
simulation period; sec Figure [T] (b). 

5.1.3 Long time behavior of the method 

It is well known that time-domain BEM analysis of- 
ten suffers from either strong numerical damping or 
instabilities when the calculation time period is long. 
Here, the long-time behavior of the present frequency 
domain method is studied. The simulation conditions 
are the same as the case 3 in Section 15.1.11 but the re- 
sponse period and the number of samplings are set to 
be T — 0.035s and N u = 512, respectively. Thus, the 
corresponding time step is At = 6.84 • 10~ 6 s and 257 
frequency domain BEM analyses are needed. 

The computed traction history is shown in Figure [5] 
Clearly, the proposed method shows a very nice behav- 
ior in the entire simulation period. Computed traction 
responses of similar problems can be found in Figure 
12 of [55] and Figure 9 of [SJ - It seems that the total 
numbers of time steps required in their work are much 
larger than that in the present simulation. For exam- 
ple, the number of time steps for obtaining Figure 9 of 
[5] is about 600/0.125 = 4800, whereas in the present 
method only 257 frequency steps corresponding to 257 
time steps are needed. 

5.2 Example 2: a plate with a hole subject to an 
impact loading 

The performance of the present method is further demon- 
strated by a more realistic problem; i.e., an aluminum 
plate with a hole subject to a vertical impact force as 
illustrated in Figure [6j The material properties of the 
aluminum plate are: E = 69 GPa, p = 2.7 • 10 3 kg/m 3 
and v = 0.3. For more details of the simulation condi- 
tions and the load history, the readers are referred to 
Section 4.2 of [UJ. In the present computation, a mesh 
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Fig. 2 Total number of iterations in solving the linear systems to the sampling frequencies 
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Fig. 3 Displacement response at the free end of the rod 



with 14072 triangular elements is used, T = 0.003s, 
= 160. The damping coefficient k is set to be 4, 
the SEM with K = 4 is used. The Hanning window is 
employed for frequency windowing. 




96.8 



Fig. 6 A plate with a hole (unit: mm) 



The responses of displacement and strain in the 
load-direction at points S (see Figure [S]) are computed 
using the present method and compared with the FEM 



and/or experiment results in Section 4.2 of [21]; see Fig- 
ure [7] By using the frequency windowing technique the 
Gibbs oscillations in both displacement and strain re- 
sponses in the later times are effectively removed. The 
windowed results coincide very well with the FEM re- 
sults. The agreement of the computed strain with the 
experimental measurement is also quite good overall. 

The total CPU time for this simulation is about 
5.4 hours. The memory usage is 840 MB, slightly more 
than the 723 MB in [2T|, because here we use GM- 
RES solver which is more memory consuming than the 
IDR(s) solver used in [2"T] . 



5.3 Example 3: a cube with 64 spherical cavities 

To demonstrate the performance of the presented method 
for solving large-scale elastodynamic wave propagation 
problems, a model consisting of a solid cube with 64 
identical spherical cavities was built and simulated; see 
Figure[8] The non-overlapped spherical cavities are ran- 
domly distributed inside the cube and their radius is 
0.072 m. The cube is fixed at the bottom surface (z = 
—0.5) and is subject to a uniform step traction p(t) at 
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Fig. 4 Traction response at the fixed end of the rod 
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Fig. 5 A long time traction response at the fixed end of the rod 



the upper surface (z = 0.5). The load p(t) is the same 
as that in Section I5.ll All other surfaces are traction 
free. The material is aluminum and the properties are 
the same as those defined in Section I5.2I This exam- 
ple was also used in (21 j . but here a finer mesh with 
N c = 231422 elements is used; hence, the total degrees 
of freedoms is nearly 0.7 million. The response period 
and the sampling size are set to be T = 0.008s and 
N u = 160, respectively. Damping coefficient n = 4.0 



and the SEM with K = 4 are used. The Hanning win- 
dow is employed in frequency windowing. 

The simulation takes about 101 hours, and con- 
sumes 21 GB memory. The CPU time still seems to 
be quite large. Further reduction of computational ex- 
pense can be achieved by employing a more effective 
pre-conditioning method and by parallel computing. 
The time responses of the displacement and traction 
in z-direction at four boundary points are plotted in 
Figure |H1 Compared with the ID solutions of the solid 
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Fig. 7 The displacement and strain responses at point S 




(a) Displacement 



rod presented in Section I5.ll the oscillation patterns of 
the displacement and traction seem reasonable. 



6 Summary 

Fast solvers are indispensable for large-scale wave prop- 
agation analysis. Recently, a fast frequency domain pFFT 
BEM incorporated with the exponential window tech- 
nique was developed for solving clastodynamic wave 



propagation problems [3T] . It has been shown that this 
method is stable and general in the sense that: (1) the 
method can be used to obtain accurate time responses 
for an arbitrarily chosen time period; (2) the method 
is applicable to linear elastodynamic systems regardless 
of the damping property. 

In this work, the computational efficiency of the 
method in |21| is further improved via three approaches: 
(1) Introducing large damping via the exponent window 

function. It can be proved that the damping actu- 
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Fig. 8 A unit cube containing 64 spherical cavities subject- 
ing to a step load 

ally plays the role of preconditioner, denoted as the 
artificial damping preconditioner (ADP) in this pa- 
per. Applying a large damping can in general im- 
prove the conditioning of the BEM matrices, thus 
reducing the total number of iterations required to 
achieve convergence. Numerical results show that 
the reduction is around 25% for k = 4 compared 
with k = 2. 

(2) Using frequency windowing technique in Fourier in- 
version. Exponential function with a large damping 
coefficient would cause unacceptable oscillations in 
the later time responses. The oscillations stem from 
the Gibbs phenomenon which is inherent in expand- 
ing non-periodic functions using Fourier series. The 
frequency domain windowing technique is an effec- 
tive way to abate the Gibbs phenomenon and is em- 
ployed in the Fourier inversion of the frequency do- 
main BEM results. It is found that the well-known 
Blackman and Hanning windows perform equally 
well in removing the Gibbs oscillations in the BEM 
solutions. The performance of the Hanning window 



is demonstrated by numerical results in Section [5] 
The time-domain damping in conjunction with the 
frequency-domain windowing techniques results in 
a modified Fourier transform method for elastody- 
namic wave propagation simulation. 
(3) Using solution extrapolation method (SEM) for ob- 
taining good initial guess. In this method, the initial 
guess for the current linear system is obtained from 
the solutions of the linear systems corresponding to 
several previous sampling frequencies via extrapola- 
tion. The aim is to reduce the number of iterations 
for iterative solvers. Numerical results indicate that 
the reduction can be up to 15% and is generally 
problem dependent. 

The effectiveness of the above strategies is demon- 
strated by three examples. The largest model has nearly 
0.7 million DOFs, and is simulated on a workstation 
with a Xeon 3.0GHz CPU and 30GB RAM. 
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